function [f,x,xp,y,yp,symparams,Phi] = DSGEmodel_yieldCurve(unitFree)
% This function reports the equations for the New Keynesian model with Calvo prices
% in f, and it reports the state vector (x and xp) and the control vector (y and yp)

%Define parameters
syms BETTA B CHI CHI0 THETA DELTA ALFA PHI PHIzero ZI ETA RHOR PHIpai PHIy PHIc PHIl NU U0 U0d...
     OMEGAd ...
     A_zz  A_nn     A_dd     A_pp    A_aa ...
     Kss OUTPUTss PAIss MUZss Css AA lss;
Rss = sym('Rss');
symparams=[BETTA,B,CHI,CHI0,THETA,DELTA,ALFA,PHI,PHIzero,ZI,ETA,...
     RHOR,PHIpai,PHIy,PHIc,PHIl,NU,U0,U0d,...
     OMEGAd,A_zz,A_nn,A_dd,A_pp,A_aa, ...
     Kss OUTPUTss PAIss MUZss Css AA lss,Rss];

%Define variables 
syms c_cu   r_cu   pai_cu   w_cu    l_cu     dc_cu  evf_cu   vf_cu    ...
     c_cup  r_cup  pai_cup  w_cup   l_cup    dc_cup evf_cup  vf_cup  ...
     d_ba1   c_ba1  muz_cu   epsd_cu  n_cu  paistar_cu  a_cu...
     d_ba1p  c_ba1p muz_cup  epsd_cup n_cup paistar_cup a_cup;

% Bond prices
syms       p1_cu   p2_cu   p3_cu   p4_cu   p5_cu   p6_cu   p7_cu   p8_cu   p9_cu   p10_cu ...
           p11_cu  p12_cu  p13_cu  p14_cu  p15_cu  p16_cu  p17_cu  p18_cu  p19_cu  p20_cu ...
           p21_cu  p22_cu  p23_cu  p24_cu  p25_cu  p26_cu  p27_cu  p28_cu  p29_cu  p30_cu ...
           p31_cu  p32_cu  p33_cu  p34_cu  p35_cu  p36_cu  p37_cu  p38_cu  p39_cu  p40_cu;
syms       p1_cup  p2_cup  p3_cup  p4_cup  p5_cup  p6_cup  p7_cup  p8_cup  p9_cup  p10_cup ...
           p11_cup p12_cup p13_cup p14_cup p15_cup p16_cup p17_cup p18_cup p19_cup p20_cup ...
           p21_cup p22_cup p23_cup p24_cup p25_cup p26_cup p27_cup p28_cup p29_cup p30_cup ...
           p31_cup p32_cup p33_cup p34_cup p35_cup p36_cup p37_cup p38_cup p39_cup p40_cup ...    

 
% Some static relations which we substitute into f
d_cu   = d_ba1p;

% LEVEL SPECIFICATION
%d_cu  =1+ A_dd*(d_cu-1) + d_cu^OMEGAd*epsd_cup;

% LOG SPECIFICATION
%log(d_cu) =A_dd*log(d_cu) + (OMEGAd*log(d_cu)+1)*epsd_cup);
%d_cu = exp(A_dd*log(d_cu) + (OMEGAd*log(d_cu)+1)*epsd_cup));
d_cup  =exp(A_dd*log(d_cu) + (OMEGAd*log(d_cu)+1)*epsd_cup);

% EQ 4: Firms FOC for labour
mc_cu  = w_cu /((1-THETA)*a_cu *Kss^THETA*l_cu ^-THETA);
%mc_cup = w_cup/((1-THETA)*a_cup*Kss^THETA*l_cup^-THETA);

% EQ 8: The aggregated production function to get output. 
output_cu  = a_cu *Kss^THETA*l_cu ^(1-THETA);    
output_cup = a_cup*Kss^THETA*l_cup^(1-THETA);    

% The real stochastic discount factor
mReal_cup = BETTA*d_cup/d_cu*...
    (AA*evf_cu^(1/(1-ALFA))/(vf_cup*muz_cup^((1-CHI)*(1-CHI0)) ))^ALFA...   
    *(c_cup-B*c_cu/muz_cup)^(-CHI)/((c_cu-B*c_ba1/muz_cu)^(-CHI))*muz_cup^(-CHI*(1-CHI0)-CHI0);  

%% The equations in f 
% EQ 2: Household's First-order condition for labour
if unitFree == 1
    f1= -1 + (c_cu -B*c_ba1/muz_cu)^(-CHI)*w_cu/(n_cu*PHIzero*(1-l_cu)^(-1/PHI));    
else
    f1= -n_cu*PHIzero*(1-l_cu)^(-1/PHI) + ((c_cu -B*c_ba1/muz_cu)/Css^CHI0)^(-CHI)*w_cu/Css^CHI0;
end

% EQ 3: Euler-equation for the one-period interest rate
if unitFree == 1
    f2 = -1 + mReal_cup/pai_cup*r_cu;
else
    f2 = -1 + mReal_cup/pai_cup*r_cu;
end

% EQ 5: First-order condition for prices
if unitFree == 1
    f3 = ((1-ETA)+ZI*mReal_cup*(pai_cup/PAIss^NU-1)*pai_cup/PAIss^NU*output_cup/output_cu*muz_cup ...
        - ZI*(pai_cu/PAIss^NU-1)*pai_cu/PAIss^NU)/(ETA*mc_cu) + 1;
else
    f3 = (1-ETA)+ZI*mReal_cup*(pai_cup/PAIss^NU-1)*pai_cup/PAIss^NU*output_cup/output_cu*muz_cup ...
        + ETA*mc_cu - ZI*(pai_cu/PAIss^NU-1)*pai_cu/PAIss^NU;
end

% EQ 6: The Taylor rule
%f4 = -log(r_cu/Rss) + RHOR*log(r_ba1/Rss) + (1-RHOR)*(PHIpai*(log(pai_cu)-(log(PAIss) + log(paistar_cu)) ) ...
%         + PHIy*log(output_cu/OUTPUTss)...
%         + PHIc*(dc_cu-log(MUZss)) ...
%         + PHIl*(log(l_cu/lss)) );
f4 = -log(r_cu/Rss) + (PHIpai*(log(pai_cu)-(log(PAIss) + log(paistar_cu)) ) ...
         + PHIy*log(output_cu/OUTPUTss)...
         + PHIc*(dc_cu-log(MUZss)) ...
         + PHIl*(log(l_cu/lss)));

% EQ 7: Output
if unitFree == 1
    f5 = - output_cu*(1-ZI/2*(pai_cu /PAIss^NU-1)^2)/(c_cu  + DELTA*Kss) + 1;
else
    f5 = - output_cu*(1-ZI/2*(pai_cu /PAIss^NU-1)^2) + (c_cu  + DELTA*Kss);
end

% Consumption growth
f6 = -dc_cu + log(c_cu/c_ba1)+log(muz_cu);

% The expected value of the value function 
if unitFree == 1
    f7 = -1 + (vf_cup*muz_cup^((1-CHI)*(1-CHI0))/AA)^(1-ALFA)/evf_cu;
else
    f7 = -evf_cu + (vf_cup*muz_cup^((1-CHI)*(1-CHI0))/AA)^(1-ALFA);
end

% The expression for the value function (minus)
if unitFree == 1
    f8 = -1 - (U0+d_cu*(1/(1-CHI)*((c_cu -B*c_ba1/muz_cu)/Css^CHI0)^(1-CHI)+U0d + n_cu*PHIzero/(1-1/PHI)*(1-l_cu )^(1-1/PHI)) )/vf_cu ...
        + BETTA*AA*evf_cu ^(1/(1-ALFA))/vf_cu;
else
    f8 = -vf_cu - (U0 + d_cu*(1/(1-CHI)*((c_cu -B*c_ba1/muz_cu)/Css^CHI0)^(1-CHI)+U0d + n_cu*PHIzero/(1-1/PHI)*(1-l_cu )^(1-1/PHI)) ) ...
        + BETTA*AA*evf_cu ^(1/(1-ALFA));
end

% Processes with time-varying vol
%The preference shock
f9 = -log(d_ba1p) + A_dd*log(d_ba1) + (OMEGAd*log(d_ba1)+1)*epsd_cu;


%% Nominal Bond Prices 
f10= -1 + (1/r_cu)/p1_cu;
f11= -1 + mReal_cup/pai_cup*p1_cup/p2_cu;
f12= -1 + mReal_cup/pai_cup*p2_cup/p3_cu;
f13= -1 + mReal_cup/pai_cup*p3_cup/p4_cu;
f14= -1 + mReal_cup/pai_cup*p4_cup/p5_cu;
f15= -1 + mReal_cup/pai_cup*p5_cup/p6_cu;
f16= -1 + mReal_cup/pai_cup*p6_cup/p7_cu;
f17= -1 + mReal_cup/pai_cup*p7_cup/p8_cu;
f18= -1 + mReal_cup/pai_cup*p8_cup/p9_cu;
f19= -1 + mReal_cup/pai_cup*p9_cup/p10_cu;
f20= -1 + mReal_cup/pai_cup*p10_cup/p11_cu;
f21= -1 + mReal_cup/pai_cup*p11_cup/p12_cu;
f22= -1 + mReal_cup/pai_cup*p12_cup/p13_cu;
f23= -1 + mReal_cup/pai_cup*p13_cup/p14_cu;
f24= -1 + mReal_cup/pai_cup*p14_cup/p15_cu;
f25= -1 + mReal_cup/pai_cup*p15_cup/p16_cu;
f26= -1 + mReal_cup/pai_cup*p16_cup/p17_cu;
f27= -1 + mReal_cup/pai_cup*p17_cup/p18_cu;
f28= -1 + mReal_cup/pai_cup*p18_cup/p19_cu;
f29= -1 + mReal_cup/pai_cup*p19_cup/p20_cu;
f30= -1 + mReal_cup/pai_cup*p20_cup/p21_cu;
f31= -1 + mReal_cup/pai_cup*p21_cup/p22_cu;
f32= -1 + mReal_cup/pai_cup*p22_cup/p23_cu;
f33= -1 + mReal_cup/pai_cup*p23_cup/p24_cu;
f34= -1 + mReal_cup/pai_cup*p24_cup/p25_cu;
f35= -1 + mReal_cup/pai_cup*p25_cup/p26_cu;
f36= -1 + mReal_cup/pai_cup*p26_cup/p27_cu;
f37= -1 + mReal_cup/pai_cup*p27_cup/p28_cu;
f38= -1 + mReal_cup/pai_cup*p28_cup/p29_cu;
f39= -1 + mReal_cup/pai_cup*p29_cup/p30_cu;
f40= -1 + mReal_cup/pai_cup*p30_cup/p31_cu;
f41= -1 + mReal_cup/pai_cup*p31_cup/p32_cu;
f42= -1 + mReal_cup/pai_cup*p32_cup/p33_cu;
f43= -1 + mReal_cup/pai_cup*p33_cup/p34_cu;
f44= -1 + mReal_cup/pai_cup*p34_cup/p35_cu;
f45= -1 + mReal_cup/pai_cup*p35_cup/p36_cu;
f46= -1 + mReal_cup/pai_cup*p36_cup/p37_cu;
f47= -1 + mReal_cup/pai_cup*p37_cup/p38_cu;
f48= -1 + mReal_cup/pai_cup*p38_cup/p39_cu;
f49= -1 + mReal_cup/pai_cup*p39_cup/p40_cu;

% % The one-period Real Bond price
% if unitFree == 1
%     f50 = - 1 + mReal_cup/p1Real_cu;    
% else
%     f50 = - p1Real_cu + mReal_cup;    
% end

%% The links
% For lagged interest rates
% if unitFree == 1
%     f51 = -1 + r_cu/r_ba1p;
% else
%     f51= -r_ba1p + r_cu;
% end

% For lagged consumption
if unitFree == 1
    f50 = -1 + c_cu/c_ba1p;
else
    f50 = -c_ba1p + c_cu;
end

%% The linear shocks
%For non-stationary technology
f51 = -log(muz_cup/MUZss)    + A_zz*log(muz_cu/MUZss);

%The preference shock
f52 = -epsd_cup;

%For labor shocks
f53 = -log(n_cup)      + A_nn*log(n_cu);

%Inflation target
f54 = -log(paistar_cup) + A_pp*log(paistar_cu);

%stationary technology shocks
f55 = -log(a_cup)  + A_aa*log(a_cu);

%Creating function f
f =[ f1; f2; f3; f4; f5; f6; f7; f8; f9;f10;...
    f11;f12;f13;f14;f15;f16;f17;f18;f19;f20;...
    f21;f22;f23;f24;f25;f26;f27;f28;f29;f30;...
    f31;f32;f33;f34;f35;f36;f37;f38;f39;f40;...
    f41;f42;f43;f44;f45;f46;f47;f48;f49;f50;...
    f51;f52;f53;f54;f55];

% Define the variables which needs to be log-transformed
y  = [r_cu   c_cu   pai_cu   w_cu    l_cu    evf_cu   vf_cu   ...
      p1_cu   p2_cu   p3_cu   p4_cu   p5_cu   p6_cu   p7_cu   p8_cu   p9_cu   p10_cu ...
      p11_cu  p12_cu  p13_cu  p14_cu  p15_cu  p16_cu  p17_cu  p18_cu  p19_cu  p20_cu ...
      p21_cu  p22_cu  p23_cu  p24_cu  p25_cu  p26_cu  p27_cu  p28_cu  p29_cu  p30_cu ...
      p31_cu  p32_cu  p33_cu  p34_cu  p35_cu  p36_cu  p37_cu  p38_cu  p39_cu  p40_cu ];
yp = [r_cup  c_cup  pai_cup  w_cup   l_cup   evf_cup  vf_cup  ...
      p1_cup  p2_cup  p3_cup  p4_cup  p5_cup  p6_cup  p7_cup  p8_cup  p9_cup  p10_cup ...
      p11_cup p12_cup p13_cup p14_cup p15_cup p16_cup p17_cup p18_cup p19_cup p20_cup ...
      p21_cup p22_cup p23_cup p24_cup p25_cup p26_cup p27_cup p28_cup p29_cup p30_cup ...
      p31_cup p32_cup p33_cup p34_cup p35_cup p36_cup p37_cup p38_cup p39_cup p40_cup ];
x  = [d_ba1  c_ba1  muz_cu            n_cu   paistar_cu   a_cu];
xp = [d_ba1p c_ba1p muz_cup           n_cup  paistar_cup  a_cup];


% Make f a function of the logarithm of the state and control vector
f = subs(f, [x,y,xp,yp], exp([x,y,xp,yp]));

% Re-defining the controls and the states
y  = [c_cu    r_cu    pai_cu   w_cu    l_cu    dc_cu   evf_cu   vf_cu    ...
      p1_cu   p2_cu   p3_cu   p4_cu   p5_cu   p6_cu   p7_cu   p8_cu   p9_cu   p10_cu ...
      p11_cu  p12_cu  p13_cu  p14_cu  p15_cu  p16_cu  p17_cu  p18_cu  p19_cu  p20_cu ...
      p21_cu  p22_cu  p23_cu  p24_cu  p25_cu  p26_cu  p27_cu  p28_cu  p29_cu  p30_cu ...
      p31_cu  p32_cu  p33_cu  p34_cu  p35_cu  p36_cu  p37_cu  p38_cu  p39_cu  p40_cu];
yp = [c_cup  r_cup pai_cup  w_cup   l_cup   dc_cup  evf_cup  vf_cup   ...
      p1_cup  p2_cup  p3_cup  p4_cup  p5_cup  p6_cup  p7_cup  p8_cup  p9_cup  p10_cup ...
      p11_cup p12_cup p13_cup p14_cup p15_cup p16_cup p17_cup p18_cup p19_cup p20_cup ...
      p21_cup p22_cup p23_cup p24_cup p25_cup p26_cup p27_cup p28_cup p29_cup p30_cup ...
      p31_cup p32_cup p33_cup p34_cup p35_cup p36_cup p37_cup p38_cup p39_cup p40_cup ];
x  = [d_ba1   c_ba1  muz_cu  epsd_cu  n_cu   paistar_cu  a_cu];
xp = [d_ba1p  c_ba1p muz_cup epsd_cup n_cup  paistar_cup a_cup];


Phi = [];
end